X-ray fluorescence spectroscopy (XRF) for metallome analysis of herbarium specimens

Background “Herbarium X-ray Fluorescence (XRF) Ionomics” is a new quantitative approach for extracting the elemental concentrations from herbarium specimens using handheld XRF devices. These instruments are principally designed for dense sample material of infinite thickness (such as rock or soil powder), and their built-in algorithms and factory calibrations perform poorly on the thin dry plant leaves encountered in herbaria. While empirical calibrations have been used for ‘correcting’ measured XRF values post hoc, this approach has major shortcomings. As such, a universal independent data analysis pipeline permitting full control and transparency throughout the quantification process is highly desirable. Here we have developed such a pipeline based on Dynamic Analysis as implemented in the GeoPIXE package, employing a Fundamental Parameters approach requiring only a description of the measurement hardware and derivation of the sample areal density, based on a universal standard. Results The new pipeline was tested on potassium, calcium, manganese, iron, cobalt, nickel, and zinc concentrations in dry plant leaves. The Dynamic Analysis method can correct for complex X-ray interactions and performs better than both the built-in instrument algorithms and the empirical calibration approach. The new pipeline is also able to identify and quantify elements that are not detected and reported by the device built-in algorithms and provides good estimates of elemental concentrations where empirical calibrations are not straightforward. Conclusions The new pipeline for processing XRF data of herbarium specimens has a greater accuracy and is more robust than the device built-in algorithms and empirical calibrations. It also gives access to all elements detected in the XRF spectrum. The new analysis pipeline has made Herbarium XRF approach even more powerful to study the metallome of existing plant collections. Supplementary Information The online version contains supplementary material available at 10.1186/s13007-022-00958-z.


Background
Herbaria around the world are the greatest sources of taxonomic, genetic, and biogeographic information on plants. As stores of well-curated plant material, they also represent an opportunity to gather further information via non-destructive techniques, but this potential remains largely unexplored to date. Recently, we have proposed the "Herbarium X-ray Fluorescence (XRF) Ionomics" approach to quantify elemental concentrations from herbarium-based plant collections [1]. Ionomics is a term referring to the study of the complete ionomic or elemental composition of a plant species [2,3] and is congruent with the metallome or elementome (sensu [4]) in this context. It requires technologies that enable high-throughput elemental analysis on large numbers of samples to gain insights in how foliar elements correlate to the ecophysiology of different plant species and how they are regulated. The herbarium XRF ionomics approach enables to quantify most elements ranging from aluminium to uranium in less than one minute [5]. This cutting-edge methodology is game-changing for the analysis of element concentrations in plant leaves across phylogenetic lineages [1].
The XRF instrument subjects a sample with a beam of high-energy X-rays generated from a X-ray tube inside the device which causes X-ray fluorescence to occur in the sample [6]. These excited fluorescent X-rays are then recorded and analysed to calculate the relative concentrations of elements present in the sample [7]. The technique has low detection limits for many different elements of the Periodic Table [8]. Compared to wet chemical analysis that involves cumbersome acid-digestion followed by ICP-AES/MS analysis, handheld XRF is a rapid method elemental analysis of plant material, and is therefore appealing to plant scientists when analysing hundreds or even thousands of samples [9][10][11][12]. In combination with simple sample grinding, a handheld XRF instrument can achieve an accuracy comparable to digestion-based methods with inductively coupled plasma atomic emission spectroscopy (ICP-AES) and surpases the latter for refractory elements, such as silicon or chromium, because of incomplete solubilization during the digestion process [10]. Many studies have developed a single element calibrations for plant matrices [9,11], but a new calibration is still required for quantifying an element that is not included in the existing calibration dataset [13].
When performing XRF scanning of herbarium specimens, only 30 s is needed to obtain quantitative data [5], it does not cause any damage the herbarium specimens [1], and it reduces the resources spent on collecting specimens [21], thereby enabling high throughput systematic scanning of herbarium collections orginating from around the world. The herbarium XRF approach was first conducted at the Forest Research Centre (Sabah Forestry Departement) in Sabah (Malaysia) on the genus Antidesma, and a critically endangered Ni hyperaccumulator plant species was identified [22]. A follow-up study led to the discovery of 28 Ni hyperaccumulator plant species, 12 Co hyperaccumulator plant species, and 51 Mn hyperaccumulator plant species from the 7300 specimens analyzed [5]. No hyperaccumulator plant species were known from Papua New Guinea before scanning of herbarium specimens of native plants from that country identified 10 Zn hyperaccumulator plant species, 8 Mn hyperaccumulator plant species, 1 Ni hyperaccumulator plant species [23]. Furthermore, herbarium XRF scanning was performed in New Caledonia on 11,200 specimens and numerous new hyperaccumulator plant species were identified: 63 Mn hyperaccumulator plant species, 5 Co hyperaccumulator plant species, 34 Ni hyperaccumulator plant species, and 4 Zn hyperaccumulator plant species [24]. In Central America, herbarium specimens of Psychotria costivenia and P. grandis were identified to be Ni hyperaccumulator plant species using XRF scanning of herbarium specimens [25].
Three different approaches are usually employed to convert the measured intensities of fluorescence X-rays in the XRF spectrum to elemental concentrations in the sample, namely Fundamental Parameters, Compton normalization, and Empirical Calibration [26,27]. A relationship between measured fluorescence peak intensities and the concentration of elements in the sample is theoretically expressed as a complex physics equation [28]. Fundamental Parameters is an approach that solves the equation by providing all parameters required by the equation [29]. In the XRF quantification process, a sample is classified into thin sample, intermediate sample, or thick sample. For thin samples, both enhancement and absorption effects are minimized, and the analyte signals follow a linear function of thickness [30]. Thick samples exceed a saturation thickness above which the intensity of the characteristic lines is constant and independent of sample thickness [30], as illustrated by Fig. 1a, b Between these extremes, the fluorescence intensities vary with thickness ( Fig. 1c-f). The handheld XRF instrument is usually calibrated based on Fundamental Parameters and assumes the sample to meet the requirement of an infinitely thick sample [31]. Because of that, it requires samples to have a certain thickness and to be homogenized (i.e., by grinding, mixing, and compacting) [21]. Failing to fulfil these requirements will cause systematic errors [26]. The Compton normalization technique normalizes the intensities of fluorescence peaks to Fig. 1 The effects of thickness on the XRF readings. Red lines are a model fitted into the reported concentrations (black dots) following Eq. (1). The samples (black dots) are constant in concentration but vary in thickness. The model reaches a plateau at a certain thickness that signifies the escape depth or critical thickness of the respective elements and shifts towards the right direction of the x-axis as the atomic number of the elements gets higher. Small variations in the sulphur and potassium XRF results are due to random noise and indicate that the samples are considered as thick samples for XRF. The underlying data used is provided in Additional file 1 the intensity of the Compton scatter [32]. However, Compton normalization may not correct all matrix effects [33].
In the case of herbarium XRF scanning, empirical calibrations have been employed to 'correct' the raw XRF-instrument generated values by measuring a set of dry leaves and establishing a regression model between these XRF values and ICP-AES measured values of these same samples [22][23][24][25]. Empirical calibration offers a simple and practical solution to correct the XRF result without any information on the properties of the instrument and sample required, but it is prone to errors [13] because it expects the calibration standards to have the exact same matrix as the actual samples which may vary in thickness, density, and composition [34]. Compared to the empirical calibration, Fundamental Parameters approaches, such as Dynamic Analysis implemented in GeoPIXE software, do not require standards, and can be used for predicting the concentration of any element, while the empirical method needs a set of standards for each different element and sample type [29,[35][36][37]. However, the major shortcoming of the Fundamental Parameters approach is that prior knowledge of the sample areal density (g/cm 2 ) is required [38,39]. Many approaches to simultaneously determine the areal density of a sample and its elemental composition have been proposed, including adding an internal standard [40] or combining XRF with other techniques [41], however none of these methods is compatible with herbarium XRF measurements, which must be fully non-destructive. Alternatively, the areal density of a sample can also be determined using the emission-transmission technique [42], which requires three measurements: that of a herbarium specimen, that of a substrate (e.g. a pure thick element), and that of a herbarium specimen on top of that substrate [35,42]. Dynamic Analysis (DA), which is based on Fundamental Parameters, has been developed for nuclear microprobe analysis and was extended for use in synchrotron-based XRF analysis [43,44]. It first uses a non-linear least-squares approach to fit the non-linear parameters, such as peak width, energy calibration and tailing in the spectrum, and then performs a further linear fit iteration with the non-linear parameters fixed in order to construct a transform matrix (the DA matrix). In contrast to early Fundamental Parameter approaches [45], this transform is applied directly to spectral data to deconvolute the spectrum into the fluorescence components for each element [46]. As shown in Fig. 2, use of the DA transform then acts like a linear LS fit. One of the many benefits of this approach, aside from fast execution, is that it can deal with complex spectra with numerous line overlaps, for example, the known interference between the Co Kα 1 line (at 6.9303 keV) and the Fe Kβ 1 line (at 7.058 keV) which makes it difficult to resolve Co if Fe is present at high concentrations. Figure 2 shows an XRF spectrum of the Ni-Co hyperaccumulator Glochidion sericeum, and as can be seen from Fig. 2b, Dynamic Analysis can separate the line overlaps successfully. As such, Dynamic Analysis is a powerful alternative to the built-in algorithms of a handheld XRF instrument while providing complete control over the input parameters, as opposed to the 'blackbox' proprietary software used by XRF instrument manufacturers. GeoPIXE is a software package that implements Dynamic Analysis, and provided that all of the required information for the Fundamental Parameters can be obtained, the calculated sample elemental concentrations are fully quantitative [47].
This article presents and validates a new data analysis pipeline (based on Dynamic Analysis) to process the raw data obtained from handheld XRF instruments for the measurement of the ionome/metallome of herbarium specimens. This pipeline is aimed to be universal permitting full control and transparency throughout the quantification process independent of the XRF instrument make or model. Furthermore, we will combine this pipeline with a new method based on the emissiontransmission technique to determine the areal density of herbarium specimens. The new pipeline will be compared to the results obtained from the XRF instrument built-in algorithms and to the empirical calibration approach.

Results
The Dynamic Analysis approach requires a complete description of both the instrument and the sample. However, many properties of the instrument are proprietary and are often undisclosed by the manufacturers. However, even with incomplete knowledge of the XRF system, it is possible to refine the model using simple standards. In this case, three 100 nm thin films [Ti, gold (Au), and tin (Sn)] with well-spaced lines covering multiple energy regions, were used to refine the instrument's sensitivity and absorption characteristics. The thin films were measured by the handheld XRF device used in this study, and the raw spectra of the thin films were processed in GeoPIXE. Parameters related to the instrument, such as the X-ray source, distance, angle, filter, and detector were based on available information as well as physical inspection of the device, and then refined using the thinfilm standard measurement data to yield results close to the expected 100% concentration for each element. Table 1 shows the final concentrations of Ti, Au, and Sn obtained in this study, and the used parameters can be found in Additional file 1.
Herbarium specimens are typically affixed to a cardboard backing, and this must be dealt with during the measurement. The composition and thickness of this cardboard may vary from herbarium to herbarium. An example of such a cardboard backing was measured while placing a 2 mm thick pure Ti plate underneath. Artefact peaks were observed between 6 and 13 keV as shown in Fig. 3c. ICP-AES analysis confirms that the cardboard contains Ca at 110, 115 µgg −1 , and this is the likely source of the strong Ca peak. The ICP-AES results also confirm that the cardboard contains Fe, but only at 69 µg g −1 , which is below the stated detection limit of the handheld device. To investigate the origin of these peaks, the 'pure' excitation spectrum of the device was obtained by triggering the X-ray tube of the handheld device into the SDD detector of a benchtop micro XRF instrument (IXRF ATLAS X). This confirmed the composition of the anode (Ag) and showed that the artefact peaks were not present in the incident beam itself, but most likely originate from secondary fluorescence from elements (such as metal alloy components) in the chamber of the handheld spectrometer and/or diffraction effects (which cannot be determined easily).
To generate a training and testing dataset, a total of leaves from 588 'normal' plants and hyperaccumulator plants were cut into 6 mm diameter discs and measured using the handheld XRF instrument. Empirical models of the areal density were established based on Ti fluorescence radiation originating from the Ti plate under the sample. Two models were generated. The first model is shown in Fig. 4 (black line) and generated by fitting an exponential model into a scatterplot between the Ti Fig. 2 Fitted XRF spectrum of a sample processed in GeoPIXE in a non-linear least-square fit (black trace) to an XRF spectrum of Glochidion sericeum (red trace) in a. Glochidion sericeum is a Ni-Co hyperaccumulator plant. The fitted spectrum is deconvolved into single fluorescence element peaks b. Note that Dynamic Analysis can deconvolve overlapping peaks such as the Fe Kβ peak from the Co Kα peak, and the Mn Kβ peak from the Fe Kα peak concentrations produced by giving the areal density values of 10.56 mg/cm 2 (mean areal density of all training dataset) to all XRF spectra and the apparent Ti concentrations produced using the measured areal density for each corresponding sample. The measured areal density was obtained by dividing the mass of each sample with the size of 6 mm diameter circle. The second model, as shown in Fig. 4 (blue line), is produced by fitting an empirical model into a scatterplot between the apparent Ti concentration produced by using the measured areal density of each corresponding against the values of the measured areal density. By using the two models subsequently, the mean absolute error of the predicted areal density is 3.19 mg/cm 2 for training dataset and 3.97 mg/ cm 2 for the testing dataset. Figure 5 shows the flowchart of the procedure used in this study to quantify elemental concentrations in GeoPIXE. The composition of (herbarium) plant leaves was assumed to be that of cellulose (C 6 H 10 O 5 ). Note that the herbarium cardboard is considered to be part of the the actual plant leaf for the purpose of the XRF analysis, but given that the cardboard is usually 'clean' for all elements of interest this has no further relevance [13]. These training samples were also used to create empirical calibration models. The 6 mm rounds were analysed by ICP-AES to determine absolute elemental concentrations. Then, a regression model between the ICP-AES results and the results of the instrument's built-in algorithm was developed. Figure 6 shows elemental concentrations from the training dataset, with the ICP-AES  Models for predicting the areal density of unknown samples, using the degree of attenuation of the signal from a pure Ti plate underneath each sample. First, a fixed areal density is used to produce an initial Ti concentration (black line). This Ti concentration is then used to predict the areal density of the specific sample via (blue line). The initial areal density used to produce dataset on the x-axis of (black line) is 10.56 mg/cm 2 which is the mean measured areal density of the training dataset results plotted against both the empirical calibration and the Dynamic Analysis. Kruskal-Wallis tests were performed to assess whether the predicted results are similar or not to ICP-AES results, and the results indicated that only the concentrations of K and Co predicted by empirical calibration are statistically similar to the ICP-AES results. The performance of Dynamic Analysis was also compared to the built-in algorithm and empirical calibration in predicting the elemental concentration of the testing dataset as shown in Fig. 7. Empirical calibration and Dynamic Analysis performed better than the built-in algorithm with errors up to 15% of the built-in algorithm errors in predicting potassium (K), calcium (Ca), iron (Fe), Co, and Zn concentrations. However, errors generated by empirical calibration and Dynamic Analysis methods were comparable or worse with 30% more errors than the built-in algorithm in predicting the concentrations of Mn and Ni (Fig. 7b). Also, from Fig. 7c it can be seen that the empirical calibration generates negative values for the concentrations of Ca, Fe, and Ni (consistent with the results based on the training dataset, see Fig. 6h).
Unlike the empirical calibration that requires a calibration line for each element, the Fundamental Parameter approach, as implemented in Dynamic Analysis, uses a global calibration which can then be applied to all elements, provided the lines used for the calibration are well-spaced. To demonstrate this, we used our built model to reprocess a total of 30,013 raw XRF spectra obtained in previous studies in Sabah (Malaysia) (7377 specimens) [23], New Caledonia (11,200 specimens) [24], and Queensland (11,436 specimens). The Dynamic Analysis approach identified several herbarium specimens with high concentrations of selenium, arsenic, and yttrium, which were not identified by the built-in algorithm and the empirical calibration (Fig. 8) and provides a good estimate of their relative concentrations in the specimens.
The heterogeneity within the specimens reduces the accuracy of the XRF measurement, and averaging multiple measurements is recommended to increase the accuracy [48]. Twenty three intact leaves were analyzed with µXRF analysis to produce high-resolution elemental maps of the leaves showing spatial variation in elemental concentrations. The average elemental concentrations of the whole leaves were compared to the concentrations of subset areas to assess the effect of the number of measurements and of the location on accuracy. In each leaf, three subset areas with a shape of a 6 mm-diameter circle (matching the XRF measurement area of the leaf discs used earlier) were drawn around the apex (1), lamina (2), and base (3) as shown in Fig. 9. Mean absolute percentage error concentrations between the whole leaf and the subsets or their combinations are calculated and tabulated in Table 2. The average of three measurement at the apex, lamina, and base consistently attains the best accuracy ~ 5% or lower with an average of 2.3% across eight different elements (K, Ca, Mn, Fe, Co, Ni, Zn, and Se). Only one measurement at apex or base generates Red hatches indicate non physical/negative values errors more than 20% for K and around 9% on average. Meanwhile, one measurement at lamina performs twice better than one measurement at the apex or base (4.3% compared to 9.5% and 9.1%, respectively). This number is comparable or better than two combinations of the three (apex:lamina:1-2 at 4.2%, apex-base:1-3 at 3.7%, laminabase:2-3 at 5.6%).

Discussion
Three 100 nm films made of pure Ti, Au, and Sn metals deposited on Kapton film were used to deduce the instrument-specific parameters, related to the detector (diameter, distance, thickness, size, and tilt angle), filters (material and thickness), and source (filter material and thickness). It was achieved by adjusting the parameters to yield Ti, Au, and Sn concentration as close as possible to 100%. After trial and error, the optimum results are shown in Table 1 with errors less than 5% for all fluorescence lines. During inspection of the raw spectra, artefact peaks were observed in the vicinity of Ti K line, Au M and L lines, and Sn L line peaks (Fig. 3b). To investigate whether the origin of these artefact peaks are from the built-in X-ray source, the X-rays emitted by the X-ray source were measured with an external XRF detector without any filter between the X-ray source and the detector except for the air path. The results show that there are no artefact peaks in the obtained spectrum, thus eliminating the possibility of the artefact peaks arising from the X-ray source. Simple background subtraction was performed to remove these artefact peaks. However, the Dynamic Analysis algorithm did not recognize the subtracted spectra because several channels or bins of the subtracted spectrum have a negative value. Due to the complexity of this issue, no functionality to remove the artefact peaks is available yet in GeoPIXE. However, it may be possible to implement a facility to make a fake "element" with peak energies, widths and relative intensities that can be specified to address this issue. In the case of spectrum range in the vicinity of Sn K line peaks, no artefact peaks were present, and the reported concentration derived from Sn K line peaks was 100%. As Fundamental Parameters requires all parameters included in the calculation, failing to correct the artefact peaks may hamper the Dynamic Analysis approach in obtaining 100% accuracy.
The built-in algorithm performs best at predicting the concentrations of Mn and Ni with the errors of the empirical calibration and Dynamic Analysis up to 24% more than the built-in algorithm errors (Fig. 7). Considering the positions of Mn and Ni fluorescence peaks, the artefact peaks from the instrument may be the reason for poor performance of both the empirical calibration and the built-in algorithm as they are in the same vicinity. For other elements, the higher errors produced by the built-in algorithm are not surprising: the manufacturer algorithms are designed for hard/dense materials (such as rocks and metal alloys) which are very different in composition, density and thickness compared to herbarium specimens [13]. The empirical calibration produces negative values for Ca, Fe, and Ni concentrations, making it unreliable prediction method for this study. These negative values could be due to errors inherited from the built-in algorithm because the concentrations reported by the instrument are not corrected for the matrix effects (absorption and enhancement) and thickness variations. Furthermore, the Dynamic Analysis approach shows its capability to compensate for the absorption and enhancement effects as shown in Fig. 6f, g where the predicted concentrations are closer to the regression line. This result signifies that the proposed approach in predicting areal density of unknown samples performs as well as Dynamic Analysis with the calculated areal density.
The leaf discs used in this study are similar to intact herbarium plant leaves and as a result, sample heterogeneity and thickness affect the accuracy of the measurement. The use of the emission-transmission technique helps to slightly improve the R 2 values (Fig. 6), but the effects of sample heterogeneity (elemental distributions and thickness variations) are still not fully corrected, and this affects especially light elements (K and Ca) more because of their vertical elemental distributions within the leaf. The leaf discs prepared for this study consist of not only lamina but also veins and midribs. The thick leaf structures (veins and midribs) contribute more signals than the thin leaf structure (lamina), as illustrated in Fig. 1. It is also possible that an element is more concentrated in the thick leaf structure, thereby producing a stronger signal. A sheet of cardboard paper was placed underneath the leaf discs for the measurement, mimicking the way herbarium specimens are mounted on cardboard sheets. The ICP-AES results show that the paper contains appreciable Ca and Fe. The leaf disc may absorb most of this Ca and Fe fluorescence coming from the cardboard, but a portion will egress through the leaf disc as a dry leaf is not infinitely thick for Ca and Fe [13]. This can be seen in Fig. 6b, d. The ICP-AES analysis reports low values but, the built-in and Dynamic Analysis algorithms report high values for both elements because the algorithms quantifies the signal from the cardboard.  The throughput of the XRF scanning of herbarium specimens is a trade-off between speed, accuracy and sensivity. When targeting to obtain datasets of tens of thousands of herbarium specimens, obtaining an error of less than 5% is acceptable. Based on Table 2, overall errors less than 5% can be achieved by measuring once at the lamina (4.3%), an average of two measurements at the apex and lamina (4.2%), or at the apex and base (3.5%) or an average of three measurements at the apex, base, and lamina (2.3%). Measuring twice will slightly reduce the errors up to 0.8% compared to measuring only once at lamina, and may be practical when the number of specimens are not that great (e.g., less 1000 specimens). Meanwhile, three measurements has errors at 2.3%, but takes three times longer. Considering that the herbarium ionomics approach aims to scan millions of specimens around the world, sacrificing 2% accuracy from 2.3 to 4.3% is in our view acceptable in return to reduce the time of the measurement by 67%.

Conclusions
Herbarium XRF Ionomics is leading the way for increased discovery of hyperaccumulator plant species and the emerging field of studying the plant metallome/ elementome. The Herbarium XRF scanning approach is another value-adding proposition for maintaining collections in global herbaria [1]. We developed a new method to determine the areal density of a dry leaf by combining the emission-transmission method with Dynamic Analysis, and by using this approach, this new pipeline outperforms the built-in algorithms. This pipeline is universal and transferable to other XRF istruments from various manufacturers, provided the relevant instrument-specific parameters are known. The XRF instrument built-in algorithms are not customizable and are unable to report some elements. The GeoPIXE-based pipeline focused on the first-row transition metal and in future may be tested on solving the overlapping fluorescence lines of high Z elements with low Z elements, thus helping in discovering more types of hyperaccumulators plants. This is particularly relevant for rare earth elements, as the XRF instrument cannot excite the Kα-lines of the REEs because their absorption edges too high, and only the L-lines will hence be excited. However, the L-lines of the REEs range from 4.64 to 8.71 keV and overlap with the Kα-lines of the first-row transition metals (4.51-8.63 keV), which makes spectral fitting very challenging. Currently, GeoPIXE caters mainly for synchrotron and desktop XRF users, but given an increase in the use of handheld XRF instrumentation, further development of GeoPIXE is highly recommended to facilitate the specific requirements for handheld XRF users, for example, to address artefact line corrections and data handling.
Thus far, the potential to unlock foliar elemental information from herbarium collections has not been fully realised. Apart from focussing on trace element (hyper) accumulation, the full metallome profile could be used to infer plant origin and adaptations of traits measurable with XRF. For example, to probe the incidence of salt tolerance/halophytes on constructed phylogenies (using rubidium and strontium as proxies), to determine the incidence of carboxylate-producing cluster-roots (using yttrium as a proxy), or to study the evolution of selenium (hyper)accumulation in the giant Astragalus genus. The use a universal independent data analysis pipeline that offers full control the quantification process, as proposed in this study, will enable future investigations along these lines of inquiry.

Handheld X-ray fluorescence spectrometer
The handheld XRF insutrument (Thermo Fisher Scientific Niton XL3t 950 GOLDD+) has a built-in miniaturized X-ray Ag tube operated at 6-50 kV and with 0-200 µA max current. The instrument was set to the pre-configuration 'Soils Mode' with the 'Main filter' for all measurements. This mode incorporates Compton Normalisation to convert the measured XRF intensities into elemental Table 2 Mean absolute percentage error of mean concentrations of a subset 6 mm diameter circle around the apex (1), lamina (2), and base (3) and their combinations compared to the mean concentrations of the whole leaves  [49,50] and the escape depth of Ti fluorescence is around 3500 µm in the leaf matrix [13]. The cardboard was placed in between the sample and the Ti plate to simulate an herbarium specimen, which are always adhered to a cardboard sheet as illustrated in Fig. 3.

Areal density modelling
Herbarium specimens can be classified as an 'intermediate sample' for the K-α line of elements with atomic numbers between 16 and 43 [13]. The intensities of these elements' K-α line fluorescence in an intermediate sample are dependent on sample thickness and full matrix composition and increase with thickness increases until reaching a certain thickness that gives saturated intensity, called thick sample or critical thickness [35]. The fluorescence intensity of an element at a given thickness is defined as follows: where, I i -the intensity of ith element at given thickness, I thick i -the intensity of ith element in thick sample, ρdensity sample (g cm −3 ), t-sample thickness (cm), and µ i -total mass attenuation coefficient (cm 2 g −1 ).
The relative intensity of ith element is a linear function of the concentration of the ith element based on the fundamental algorithm [51]. Because the relative intensity of ith element in an intermediate sample is thickness dependent, as illustrated in Fig. 1, it is crucial to correct the relative intensity for thickness variations.
Dynamic Analysis, the Fundamental Parameters approach, has been developed for nuclear microprobe and synchrotron-based XRF [43,44]. The algorithm constructs a matrix transform to mimic a linear leastsquare fit to an XRF spectrum, in order to deconvolved the spectrum into fluorescence peaks for each element [46] as shown in Fig. 2. The advantages of Fundamental Parameters such as Dynamic Analysis, compared to the empirical calibrations, are that the Fundamental Parameters does not require standards, and can be used for other elements, while the empirical method needs a set (1) of standards for each element and sample type [29,[35][36][37]. The SDD detector of the handheld XRF instrument is has about 150 eV resolution. As a consequence, overlapping peaks are indistinguishable spectroscopically, for example, a common peak overlap between Fe Kβ1 line (at 7.058 keV) with the Co Kα1 line (6.9303 keV). Nevertheless, Dynamic Analysis can separate overlapping peaks of Co and Fe shown in Fig. 2b. Therefore, the Dynamic Analysis approach is a powerful alternative to the builtin algorithm of a handheld XRF instrument while having complete control over the input parameters (as opposed to the 'blackbox' propriety software used by XRF manufacturers). The GeoPIXE is a software package that is based on the Dynamic Analysis method [47].
The major shortcoming of the Fundamental Parameters approach, such as in Dynamic Analysis, is that prior knowledge of sample areal density (g/cm 2 ) is required [38,39]. Many approaches to simultaneously determine the areal density of a sample and its elemental composition have been proposed, including adding an internal standard [40] or combining XRF and other techniques [41], however none of these methods compatible with herbarium XRF scanning. The areal density of a sample can also be determined using emission-transmission technique [42]. This requires three measurements: herbarium specimen, a substrate (pure thick element), and herbarium specimen on top of the substrate [35,42]. If I s , I t , and I ts refer to the Ti intensities from the sample alone, target alone, and the sample on top of the target, respectively, then the relationship between intensities and areal density can be expressed as follows: where m , areal density (g/cm 2 ), is the product of density ρ multiplied by sample thickness t . Until now, all the herbarium XRF scanning conducting so far uses a 2 mm thick pure Ti plate [17][18][19][20], and Ti concentration in leaves is less than 34 µgg −1 or even less [49,50] (3) in areal density. Empirical models were built to quantify the variation and later could be used to predict the areal density of unknown sample. The procedure of establishing the empirical models for predicting areal density is as follows: 1. Each a 6 mm diameter leaf was placed at the geometric centre of the examination window and subsequently measured using the handheld XRF spectrometer. A stack of a sheet of cardboard, a 2 mm thick pure Ti plate, and a 2 mm thick pure Mo plate were placed below the leaf disc for the measurements. 2. The obtained XRF spectrum was processed using the Dynamic Analysis algorithm by providing an initial density 10.56 mg/cm 2 corresponding to the mean areal density of the training dataset.

The Ti concentrations reported by the Dynamic
Analysis algorithm with the initial density were plotted against Ti concentrations reported by the Dynamic Analysis algorithm with the measured areal density that were plotted also against the measured areal density.

Validation of the proposed approach
The dataset was divided into two datasets, and the proportion of the training and test dataset is 80% (470 samples) and 20% (118 samples), following the Pareto principle. The performance of Dynamic Analysis was compared against the built-in algorithm and empirical calibration approach that used the training dataset to build the empirical model. The results of ICP-AES were used as the reference data, and mean absolute error was used to quantify the magnitude of the errors generated by each approach compared to the ICP-AES results. Also, ratioing the elemental concentration of each approach to ICP-AES was performed to evaluate the direction of the errors whether underestimated or overestimated.

Dynamic analysis in GeoPIXE
Three batch processing were conducted in GeoPIXE. The first batch aimed to obtain the instrument related parameters and was achieved by processing three certified thin films (Ti, Au, and Sn). Because the parameters of the samples (the composition and thickness of the metal films) were known, the instrument related parameters in GeoPIXE were adjusted so that the concentration of the films is 100%. The second batch processing was designed Table 3 List of specimens scanned using the µXRF instrument. The letter 'X' indicates elemental concentration image available for the corresponding specimen to establish a procedure for predicting the areal density of an unknown sample. The last batch was performed on the test dataset based on the parameters and procedure obtained from the first and second batch processing. During the last batch, the continuum background of the spectrum is also corrected and estimated for fluctuated and low counts [52]. Then, the limit of detection is defined as following, 3.29 σ b , where σ b is the background standard deviation [53]. The limit of detection for each spectrum differs due to differences in matrix compositions and physical properties, and thus, the limit of detection are the mean of the detection limit. The average limits of detection for K, Ca, Mn, Fe, Co, Ni, and Zn were 2936 µg g −1 , 1829 µg g −1 , 154 µg g −1 , 127 µg g −1 , 104 µg g −1 , 134 µg g −1 , and 101 µg g −1 , respectively.

ICP-AES analysis
The weight of all samples was recorded and the samples were pre-digested in 2 mL HNO 3 (70%) for 48 h. After that, the samples were fully acid digested for 45 min in a microwave oven (Milestone Start D) using the following program rise to 80 °C within 10 min, rise to 135 °C in the next 15 min, hold for 20 min at 135 °C, and cool down for 30 min. The samples were then brought to volume (10 mL0 and analysed with an ICP-AES (Thermo Scientific iCAP 7400) instrument for natrium, magnesium, aluminium, phosphor, sulphur, K, Ca, chromium, Mn, Fe, Co, Ni, Cu, Zn, As, selenium, cadmium, and thulium in axial or radial modes depending on the expected analyte concentration. A 4-point calibration was used to calibrate all elements and to compensate for matrix-based effects, yttrium was used as an in-line internal addition standardization. Standard Reference Material (NIST Apple Leaves 1515), internal reference materials, certified reference material (Sigma-Aldrich Periodic Table mix 1 for ICP TraceCERT ® ), and matrix blanks were measured as quality controls. The results of ICP-AES analysis were used here as a reference for the XRF data calibrations.

Statistical analysis
R version 3.6.2 was used for statistical analysis. A nonlinear least square, nls, function which is the built-in function of R was used to fit an empirical model into the scatterplot between Ti concentration and areal density data. The "Metrics" package was used for calculating the mean absolute error of each prediction, and the "agricolae" package was used for conducting the nonparametric, Kruskal-Wallis, test. Meanwhile, the "ggplot" package was used for producing the charts.

Spatial variation analysis
Because herbarium specimens were measured in situ, it is important to assess the spatial variability of elemental concentrations across the leaf. In this case, twenty three leaves were scanned using a custom-built µXRF instrument at the Centre for Microscopy and Microanalysis, the University of Queensland, Brisbane, Australia to map the elemental concentrations in the whole leaves. The instrument has a 50 kV (1000 μA) Mo-target source that generates 17.4 keV X-rays fitted with polycapillary focussing optics to 25 μm. The leaf specimens were placed between two sheets of 6 μm Ultralene thin film stretched over a Perspex frame magnetically attached to the sample platform of the instrument and measured at atmospheric temperature (~ 20 °C). The results are multiband 'elemental' images and each band represents a discrete energy (and thus element). The multi-band images were processed in GeoPIXE producing elemental maps of different elements (K, Ca, Mn, Fe, Co, Ni, Zn, and Se) as shown in Fig. 9a. These elemental maps were further processed in QGIS 3.16 to segment each image into the specimen and the background (Fig. 9b). It was classified by manually setting the threshold between the specimen and the background. Three circles with a diameter of 6 mm were drawn at the apex, lamina, and base to simulate three measurements using a portable XRF instrument (Fig. 9c). The average concentration of each element per specimen was calculated (green area in Fig. 9b) and compared to the average concentration of the three spots shown in Fig. 9c. Table 3 lists specimens scanned using the µXRF instrument.